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Spherical symmetry is ubiquitous in nature. It's therefore unfortunate that simulation of spher- 
ical systems is so hard, and require complete spheres with millions of interacting particles. Here 
we introduce a method to model spherical systems, using revised periodic boundary conditions 
adapted to spherical symmetry. Method reduces computational costs by orders of magnitude, and 
is applicable for both solid and liquid membranes, provided the curvature is sufficiently small. We 
demonstrate the method by calculating the bending and Gaussian curvature moduli of single- and 
multi- layer graphene. The method works with any interaction (ab initio, classical interactions), with 
any approach (molecular dynamics, Monte Carlo), and with applications ranging from science to 
engineering, from liquid to solid membranes, from bubbles to balloons. 
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I. INTRODUCTION TO MODELING 
APPROACH 

The problem in simulating spherical symmetry is topo- 
logical: you cannot build a perfect sphere from identi- 
cal blocks. The absence of such a building block has 
enforced expensive simulations with complete spheres — 
though usually spherical simulations are simply avoided. 
Overwhelming dilemmas like this are often considered 
so fundamental and frustrating that they restrain all at- 
tempts to seek for a practical solution. 

Anyhow, avoiding spherical systems in our world is 
hard. Spherical shells surround us in a variety of forms: 
in balloons, in cell membranes inside our bodies, in bub- 
bles in the sea, or in Earth's crust. The interaction of 
nanoparticles with cell membranes, for instance, is a top- 
ical question.^ Since cell membranes' curvature moduli 
determine the very forms of red blood cells, for example, 
one can see why simulations should incorporate curvature 
effects!^ Another timely example is the foam of spheri- 
cal bubbles in the sea, the bursting of which may play an 
important role on the so-called sea sprav that produces 
spherical aerosols into the atmosphere!^ 

Although liquid membranes are more abundant in na- 
ture, also man-made solid membranes have spherical 
symmetries, at least locally. Examples are fullerenes, 6 
nanoballoonsP and especially graphene that contains in- 
trinsic ripples even when suspended freelyP^ Curvature 
moduli of graphene are intimately related to these rip- 
ples, whether they are intrinsic or not j n l 12 l and in a broad 
sense to elastic behavior of all honeycomb carbon, among 
graphene nano ribbon sP^ multilayer graphene and car- 
bon nanotubes! 14 * 15 * 

Conventionally spherical systems are treated in three 
ways. The first way is to simulate the system as a whole. 
Needless to say, this is expensive and often impossible.^ 
The second way is, should the system have some well- 
defined point-group symmetries, to use those symmetries 
for reducing computational costs. Most established codes 
have the ability to benefit from such symmetries; this 
has long been a standard procedure with molecules and 
clusters^ Because the symmetry is exact, however, nei- 
ther the curvature nor other geometrical parameters can 



be changed flexibly. The third way is to ignore curva- 
ture altogether and to use periodic boundary conditions 
(PBC) to simulate an infinitely large, flat membrane. 
Unfortunately, in nanoscience many systems fall between 
these two extremes: systems with huge number of parti- 
cles, having no overall symmetry, but prominent curva- 
ture effects. At the moment a practical way to simulate 
such systems does not exist. 

The periodic boundary conditions have been adapted, 
however, also to symmetries beyond translati on. The 
first ideas came along chiral carbon nanotubesj 17 * 18 * and 
those ideas have been used ever since; for reviews look at 
Refs. fT9landl20l An important extension to general sym- 
metries with exact treatment was done in Ref. |21j which 
has enabled more flexibility!^"^. Later, Ref. 25] intro- 
duced revised periodic boundary conditions (RPBC), a 
simple formalism for general material distortions; this is 
the approach we shall use here, and it's illustrating to 
review it briefly. 

In RPBC, the usual translation operations are replaced 
by general symmetry operations S n that, in a quantum- 
mechanical language, leave the electronic potential in- 
variant, or 

D(S n )V(r) = V(S~ n r) = V(r). (1) 

The operation S n is a succession of an abelian group 
of operations Si, that is S n = S^S^ 2 • • • • Then, by 
imposing periodicity (Sf /ti = 1, Mi integer), one finds 
that the Hamiltonian eigenstates ^ aK (r) at r and at r' = 
S~ n r differ only by a phase factor, 

D{S n )il) an {r) = ^an(S- n r) = exp(z* • n)^ a#c (r), (2) 

with inverse operation <$ _r \ band index a, and the recip- 
rocal lattice vector k. Eq.p| infers the familiar result: a 
single simulation cell — whatever its shape — is enough to 
describe the extended system as a whole. Revised PBC 
is hence similar to conventional PBC and differs only in 
the definitions of the symmetry operations. There are no 
other fundamental differences. As an illustrative exam- 
ple, the total energy with a classical pair potential is 

1 N 

^Pair = g E (3) 
i,j=l n 
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where N is the particle count and n runs over operations 
where particle i at Ri still interacts with the periodic 
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FIG. 1: (Color online) (a) Illustration of symmetry operations 
Si and S2 for spherical symmetry. «Si is a rotation of angle 
50i around y-axis and S2 is a rotation of angle 662 around 
£-axis; the angles 50i are small. (In general, c\ can also be 
non-orthogonal to C2 and 50i different from 862-) 



image of particle j at S n Rj. Forces are the negative 
gradients of this expression, as usual. Look at Ref. [25] 
for details of RPBC and Refs. [T5l and [26l for examples of 
usage. 

In this paper we use the RPBC, reviewed above, in an 
approximate way to introduce a trick for modeling spher- 
ical membranes. Adapting RPBC for spherical systems 
enables simulations with orders-of-magnitude reductions 
in computational costs. We shall apply the method to 
calculate graphene's mean and Gaussian curvature mod- 
uli, but first we proceed to discuss symmetry operations 
and their character. 



II. SPHERICITY AS AN APPROXIMATE 
SYMMETRY 

Consider the square cone in Fig. [I] regard the grid as 
fixed in space, and concentrate on the shaded region. If 
we rotate all particles an angle 50 1 around ?/-axis, or an 
angle 59 2 around x-axis, the geometry within the shaded 
region will remain approximately intact. This means that 
rotations S^r = 1Z(ni59iCi)r and S^r = 1Z(n 2 59 2 c 2 )r 
(with c\ = j, C2 = i and operation 11(c) as |c|-radian 
rotation around c) leave the electronic potential V(r) in- 
variant near the shaded region: Si and S 2 are symmetry 
operations as far as the shaded region and its vicinity is 
concerned. Two rotations around different axes do not 
commute in general, but if the rotation angles 59i are 
small Sir ~ r + 59iCi x r, rotations do commute to linear 
order in <%Vs, [<Si,c>2] = 0(59%). Hence also the com- 
bined operation 

(S^S^)r * 1Z(m50iCi + n 2 50 2 c 2 )r = S n r (4) 

is an approximate symmetry operation, provided that rti 
are small enough. Eq.Q is basically all we need to fully 



employ the RPBC of Ref. 25; we are, in principle, ready 
to go and to simulate any spherical membrane. 



III. FEATURES DUE TO APPROXIMATION 

In practice, however, the approximate character of S n 
raises questions that deserve some elaboration. First, 
as already mentioned, the formalism assumes periodic 
boundary conditions (S^ di = 1) which may seem ques- 
tionable. Here we remind that similar PBCs are used 
also in regular bulk, with all three dimensions periodic 
in an intertwined fashion. (In two dimensions PBC rep- 
resents topologically a toroid.) The bottom line is that 
periodicity is not a physical reality but a mere mathe- 
matical trick that works, and enables the application of 
revised Bloch's theorem in the first place P^ESl The inte- 
gers Mi are connected rather to K-point sampling than 
to physical reality. 

Second, revised PBC does not need the "unit cell" con- 
cept. However, we shall call the square cone in Fig. [TJ 
extending from the origin to infinity and enclosing the 
shaded region, a unit or simulation cell because the con- 
cept is familiar and convenient in discussion. Otherwise, 
the mere expression for S n in Eq.Q is enough to deter- 
mine everything in the simulation. 

Third, the claim is not to simulate a complete sphere, 
but rather to view the curvature as a local property. The 
particles in the simulation cell see the closest environ- 
ment curved — and only this is important. The simula- 
tion cell is the only cell we model, and distances and 
angles measured only from the simulation cell are mean- 
ingful. For example, the vicinity of particle at r in Fig [I] 
exhibits curvature in bond angles and distances if one 
looks at particle's own periodic images at S 2 r and S 2 ~ 1 r. 
Symmetry operations S n that have rti large enough to 
rotate large angles (rti ~ (ir/2)/59i) should be excluded 
because the non-commutativity of S^s would otherwise 
become significant. 

Fourth, the radius or curvature R in Fig. [T] is not a pa- 
rameter in the simulation; radially particles can migrate 
wherever interactions drive them. The spherical form is 
only forced by the choice of symmetry operations and 
the parameters 59 1 and 59 2l and since the symmetry is 
discrete, the system needs to be neither continuously nor 
smoothly spherical. 

Fifth, a natural limitation is to have enough empty 
space near the origin to avoid too close encounters be- 
tween the particlesP^ Membrane can be thick. 



IV. APPLYING THE METHOD 

The validity of the method depends on the system and 
its interactions. As a principal rule, the radius of curva- 
ture R should be much larger than the interaction ranges 
between the particles. If ranges are larger than the sys- 
tem size, especially if those interactions control morphol- 
ogy, one does better to model the complete system. The 
Coulomb interactions can play a role locally, within small 
length scales (size of the unit cell at most), but the long- 
ranged Coulomb interaction requires special care, per- 
haps some refinements (the unit cell better be neutral) 
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Quantitative error due to the non-commutativity of S^s 
can be estimated by first using the right-hand side of 
Eq.ffl as 5™, and then using the left-hand side of Eq.Q 
as 6^ (changing the ordering of S^S^ 2 ), and comparing 
the resulting energies. 

Because liquid lacks long-range order, the method suits 
particularly well for liquid membranes, such as lipid bi- 
layers. Their energetics can be described by the Hel- 
frich Hamiltonian that gives membrane's elastic energy 
per unit area as^ 
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Here n is the mean curvature modulus (don't confuse 
with a K-point), n is the Gaussian curvature modulus, 
and Ri and R2 are the principal radii of curvature. The 
liquid membrane doesn't need to be free-standing, how- 
ever, because also solid support can be incorporated, ei- 
ther by external force fields or by fixed atoms. External 
radial forces can be also used for pressurization, mim- 
icking the embedding of membrane in gaseous or liquid 
environments. 

For solid membranes the situation is more complicated, 
because energy will come also from the internal strain E s . 
If a flat, round sheet of radius p is wrapped into a spher- 
ical segment, the energy will be E s ~ Ehirp 6 / '108 i? 4 , 32 
where E is the Young's modulus of the material, h is the 
membrane thickness, and R is the radius of curvature; 
meanwhile the curvature-related energy is E c = g ■ 7rp 2 . 
Hence, for a reasonable modeling of solid membranes us- 
ing Eq.([5|, we need to have E s <C E c , or 



E s /E c 
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<1, 



(6) 



which suggests a minimum radius of curvature R m in for 
a given unit cell area. If this geometrical and material- 
dependent criterion should be violated, the simulation 
would be dominated by non-local stress fields. Since the 
method does not properly describe these fields, the treat- 
ment would become ill-defined. 

The above problem is present when sphericity is forced 
on originally flat sheet. But defects, for example, can in- 
duce spontaneous curvature in solid membranes in which 
case Rmi n can be smaller. The method provides a new 
tool to investigate phenomena such as rippling due to 
adsorption-induced pinching of the membraneP^ This 
method does not directly compete with any existing 
method, but instead it provides possibilities to do some- 
thing new. 



V. EXAMPLE: SPHERICAL GRAPHENE 

The spherical symmetry was implemented in the 
density- functional tight-binding code hotbit! 3 - 3 -! 3 - 4 ^ The 
RPBC implementation has a negligible computational 
overhead as compared to translational symmetry,^ and 
can be implemented just by a few lines of new code in 
any existing RPBC implementation. The code source is 
open and stands for inspection. 

In this section we use the hotbit implementation to 
present one practical example. We calculate the curva- 
ture moduli of graphene, motivated by their relevance to 
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FIG. 2: (Color online) (a) Two-atom unit cell for spherical 
graphene, illustrating the symmetry operations: Si is a ro- 
tation of angle S0i around ci and & is a rotation of angle 
5O2 around C2. (b) Few periodic images of atoms a and 6, 
shown for visualization purposes only, (c) Elastic energy per 
atom as a function of radius of curvature. Inset: fit to R~ 2 
behavior; the thin shaded fan is the error estimate due to 
approximations involved. 



present-day engineering with carbon nanostructures. For 
a sphere the radii of curvature are R\ = R2 = R, and 
Eq.([5| gives g = (2n + n)/R 2 ] for a cylinder Ri = R, 
R2 = 00, and g = k/(2R 2 ). Hence, by calculating the 
elastic energies for a cylinder and a sphere and varying 
S6i's (hence varying R) we obtain both n and n directly. 

Prior to simulating spherical graphene, we first cal- 
culated the mean curvature modulus of graphene, also 
applying revised PBC. Only now the symmetry opera- 
tions, in a cylinder-like setup, were a rotation around 
z-axis (<Si) and translation in z-direction (£2) with a 4- 
atom unit cell (like a nanotube with enormous diame- 
ter); we won't discuss the cylinder setup further hereP^ 
The resulting cohesive energy depends on R quantita- 
tively like R~ 2 , as Eq.([5| suggests, and the fitted value 
for k, = 1.61 eV (4.22 eVA 2 /atom) agrees with a density- 
functional reference value (1.5 eV)P^albeit is larger than 
an experimental reference value (1.2 eV)P^ 

Returning to spherical graphene, Figs.|2^i and|2]3 show 
the two-atom unit cell of graphene. UnliKe in Fig{l] the 
unit cell is skewed with Ci = j and C2 = cos(57r/6)£ + 
sin(57r/6)j. The geometry was optimized with given ^^'s, 
which were taken as 56i — 2.5 A/R' when we wanted 
to investigate a radius of curvature that roughly equals 
R'W* All the radii of curvature we report, anyhow, are 
the optimized R (R ~ R' because curvature changes bond 
distances only slightly). In practice we found that struc- 
ture optimizations require convergence criteria tighter 
than with translational cells, due to geometrical effects 
from small 56 iW* In quantum simulations ^-points can 
be freely sampled G [— 7r, 7r]) because PBC is an ap- 
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TABLE I: Curvature moduli for single- and multi- layer 
graphene (AB stacking). Numbers in parentheses are esti- 
mates from Eq.([7|. a ^ k — 1.610 eV for bending against 
zigzag direction (armchair direction remains straight), and 
hi — 1.606 eV for bending against armchair direction. 



layers (N) 


kn (eV) 


kn (eV) 


monolayer 


1.61 a 


-0.70 


bilayer 


180 (180) 


-140 (-176) 


trilayer 


690 (660) 


-600 (-700) 



proximation, just as with conventional Bloch's theorem; 
we used a 50 x 50 K-point mesh. 

Figj2]3 shows our main result, graphene's cohesive en- 
ergy as a function curvature — and represents the show- 
case of the new physics this method can unearth. En- 
ergy behaves clearly like R~ 2 , as suggested by Eq.([5]). 
The energy penalty 6.6 eVA 2 .R _2 /atom, combined with 
previously calculated ft, yields ~K = —0.70 eV; we could 
not find this number in the literature. This result con- 
firms graphene's beautiful elastic behavior up to high 
curvature — also for spherical distortion.^ 

We did consistency checks for the graphene sphere 
calculations, three listed next. As a first check, when 
we investigate Eq.(|6]) with graphene parameters, we get 

p <C a/6 A • R. For a graphene unit cell p ~ 1 A (lattice 
constant 2.5 A), and the consequent criterion R ^> 0.2 A 
is easily fulfilled. We obtained the same n with N = 8 
and N = 32 atom unit cells, even though larger N in- 
creases i? m in (Eq.([6| and p 2 oc N infer R m i n oc N). Thus, 
the area is small enough to be stress- free, and the sim- 
ulation is indeed dominated by curvature energy alone. 
We were able to perform controlled calculations down to 
radii R m in ~ 10 A or 50 max ~ 15°. As a second check, 
we estimated quantitative error in energy due to the non- 
commutativity of the two rotations (inset in Fig|2fc), as 
suggested above, but found the error fairly small. As 
a third check, we implemented symmetry also with a 
negative Gaussian curvature R\ = — R2 = R, for which 
g = —k/R 2 directly, and got an independent confirma- 
tion for K; we won't attempt to describe structures with 
negative Gaussian curvature here. Finally, since there is 
no charge transfer, the long-range Coulomb interactions 
are no issue. 

Closer inspection of geometry revealed that curvature 
increased bond distances as d nn = 1.417 A+0.135 A 3 /i? 2 , 
due to the weakening of in-plane a-bonds, and hereby 
decreasing the effective nearest-neighbor tight-binding 
hopping parameter as t e R = t gr — 4.8 eVA 2 /R 2 (t gr « 



2.7 eV). For a detailed discussion of curvature-induced 
effects on graphene, we recommend Refs. l39l l40l and I4T1 
For completeness we calculated n and ~R for bi- and 
trilayer graphene as well, and summarize the results in 
Table [TJ Assuming a constant layer separation of h = 
3.4 A analytical expressions for the curvature moduli of 
multi-layer graphene come as 

K n =n«i + Eh 3 (n 3 - n)/12 , 

(7) 

K n =riR\ - Eh 3 (n 3 - n)/12, 

where n is the number of layers and E is Young's mod- 
ulus. The simulated and analytical numbers have a fair 
agreement. Table reveals how strikingly smaller the mod- 
uli are for graphene monolayer, a true oddity among solid 
elastic sheets, as noted already in Ref. |42j 



VI. CONCLUDING REMARKS 

We have introduced a simple and practical method to 
simulate spherical systems using revised PBC. Although 
the method is approximate, it is applicable precisely to 
systems so hard to handle: large systems with prominent 
curvature effects. Since the method works with schemes 
from ab initio electronic structures and classical poten- 
tials to coarse-grained and finite element modeling, and 
has a wide range of applicability, we encourage any ad- 
ditional implementations. 

Admittedly, it may take some time to digest the ap- 
proximate nature of the method. The role of symme- 
tries in materials modeling is usually taken as clear-cut, 
solid, and untouchable: it either is or is not. In this pa- 
per we have, however, created and entered a new gray 
area in symmetry usage; we are unaware of symmetry 
being treated in this type of approximate fashion be- 
fore. For this reason, when using approximate spherical 
symmetry — or other approximate symmetries in future — 
we urge to examine modeled systems carefully and get 
assured of method's validity; the best guide on this way 
is common sense. 
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